
#include <msGeneralizedCoordinates.h>
#include <msParser.h>
#include <msZMat.h>
#include <msSystem.h>
#include <msKineticOperator.h>
#include <msLagrangian.h>
#include <muParser.h>
#include <msQuantumMotion1d.h>
#include <msLinearRigidRotor.h>
#include <msTranslation.h>
#include <msSolverLagrangian.h>


int main(int argc, char **argv) {
 
 using namespace atomism;

 //****************
 // Initialization
 //****************

 boost::shared_ptr<msUnitsManager> unitMng 	= 		msUnitsManager::New();

 unitMng->set( "Angstrom Degree amu hartree ps" );

 boost::shared_ptr<msZMat> zmt			=  		msZMat::New(unitMng);
 zmt->setId("zmat");
 zmt->set("H H 0 0.74191");

 boost::shared_ptr<msSystem> system 		=  		msSystem::New(unitMng);
 system->addEntity(zmt);

 boost::shared_ptr<msLinearRigidRotor> rotor 	= 		msLinearRigidRotor::New(unitMng);
 rotor->setMmtOfInertia(system);
 rotor->setSymmetryNumber(2);
 rotor->getParameters()->setParameter("Emax",200,"kcal/mol");
 
 boost::shared_ptr<msTranslation> 	trans	= 		msTranslation::New(unitMng);
 trans->setSystem(system);
 trans->setId("trans");

 boost::shared_ptr<msGeneralizedCoordinates> gcoors = 		msGeneralizedCoordinates::New(unitMng);

 boost::shared_ptr<msScalarVariable> q0 	= 		msScalarVariable::New("Angstrom");
 q0->set( 0.7414, 0.1, 20. , 0.01, 0.1);
 q0->setId("q0");
 gcoors->addVar(q0);

 double De = unitMng->convert("eV",4.7457);
 boost::shared_ptr<msParser> epot = 				msParser::New(unitMng->getEnergyStr());
 epot->setCoordinates(gcoors);
 epot->setId("f=Morse");
 epot->setExpression("De * ( 1 - exp( - beta*( q0 - re ) ))^2");
 epot->setConstant(	"De"  , De );
 epot->setConstant(	"beta", unitMng->convert("m^-1",1.94196e10) );
 epot->setConstant(	"re"  ,	unitMng->convert("m",0.74191e-10) );

 boost::shared_ptr<msParser> epotHO = msParser::New(unitMng->getEnergyStr());
 epotHO->setCoordinates(gcoors);
 epotHO->setExpression("alpha * ( q0 - re )^2");
 epotHO->setId("f=HO");
 gcoors->setCoordinate(0,	unitMng->convert("m",0.74191e-10));
 epotHO->setConstant("alpha",epot->partialSecondDerivative(0) );
 epotHO->setConstant("re",	unitMng->convert("m",0.74191e-10));

 boost::shared_ptr<msParser> kinfct = 				msParser::New("Angstrom");
 kinfct->setCoordinates(gcoors);
 kinfct->setId("f=q0");
 kinfct->setExpression("q0");

 boost::shared_ptr<msKineticOperator> kinop =  			msKineticOperator::New(unitMng);
 kinop->set(gcoors,system);
 kinop->setId("kinetic_operator");
 kinop->addDynamicDof("zmat",0,kinfct);

 boost::shared_ptr<msQuantumMotion1d> qmotion = 		msQuantumMotion1d::New(unitMng);
 qmotion->setEqOfMotion(kinop,epot);
 qmotion->getParameters()->setParameter( "BasisSize",501,"");
 qmotion->getParameters()->setParameter( "Emax",De,unitMng->getEnergyStr());
 
 boost::shared_ptr<msQuantumMotion1d> oscillsHO = 		msQuantumMotion1d::New(unitMng);
 oscillsHO->setEqOfMotion(kinop,epotHO);
 oscillsHO->getParameters()->setParameter("BasisSize",201,"");
 //*************
 // Computation
 //*************
 
 qmotion->freezeKinValue(kinop->kineticFunction());
 qmotion->computeEigenValues();
 fprintf(stdout,"Morse vibrator:\n");
 fprintf(stdout,"0th eigenvalue / exact one =%g ; \n",qmotion->eigenvalue(0)/unitMng->convert("hartree",9.86922e-3) ); 
 fprintf(stdout,"4th eigenvalue / exact one =%g ; \n",qmotion->eigenvalue(4)/unitMng->convert("hartree",7.847518e-2) ); 
 fprintf(stdout,"9th eigenvalue / exact one =%g ; \n",qmotion->eigenvalue(9)/unitMng->convert("hartree",1.3836303e-1) ); 

 fprintf(stdout,"Number of bounded states = %g ;\n", qmotion->SOS()->msScalarFunction::evaluate(unitMng->convert("eV",4.7456)) );
 
 oscillsHO->computeEigenValues(); 
 fprintf(stdout,"HO vibrator:\n");
 fprintf(stdout,"0th eigenvalue / exact one =%g ; \n",oscillsHO->eigenvalue(0)/unitMng->convert("hartree",9.86922e-3) ); 
 fprintf(stdout,"4th eigenvalue / exact one =%g ; \n",oscillsHO->eigenvalue(4)/unitMng->convert("hartree",7.847518e-2) ); 
 fprintf(stdout,"9th eigenvalue / exact one =%g ; \n",oscillsHO->eigenvalue(9)/unitMng->convert("hartree",1.3836303e-1) ); 

 qmotion->getEigenvector(0)->save(0,500,"ev0.dat");
 qmotion->getEigenvector(4)->save(0,500,"ev4.dat");
 qmotion->getEigenvector(9)->save(0,500,"ev9.dat");

 qmotion->SOS()->save(0,500,"qmotionSOS.dat");

 rotor->computeEigenValues();
 rotor->getOccupencyProbability(500)->save(0,500,"occProba_J_500.dat");
 rotor->getOccupencyProbability(2000)->save(0,500,"occProba_J_2000.dat");

 double S_300 = trans->S(300) + rotor->S(300) + qmotion->S(300) + unitMng->convert("J/mol",csts::R);
 double Cp_300 = trans->Cv(300) + rotor->Cv(300) + qmotion->Cv(300) + unitMng->convert("J/mol",csts::R);
 
 fprintf(stdout,"total with morse vibrator:\n");
 fprintf(stdout,"S(300 K) / S_janaf(300 K) = %g ; \n",S_300/ unitMng->convert("J/mol",130.858) );
 fprintf(stdout,"Cp(300 K) / Cp_janaf(300 K) = %g ; \n",Cp_300/ unitMng->convert("J/mol",28.849) );
 
 ofstream file("TD.dat");
 file<<"T_K S_"<<unitMng->getEnergyStr()<<" Cp_"<<unitMng->getEnergyStr()<<endl;
 for( double T=200;T<3000;T+=100)
    file<<T<<" "
        <<trans->S(T)  + rotor->S(T)  + qmotion->S(T)  + unitMng->convert("J/mol",csts::R)<<" "
        <<trans->Cv(T) + rotor->Cv(T) + qmotion->Cv(T) + unitMng->convert("J/mol",csts::R)<<endl;
	
 file.close();
 return 0;
 }


